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Abstract 

Simple finite differencing of the anisotropic diffusion equation, where diffusion is only along a given 
direction, does not ensure that the numerically calculated heat fluxes are in the correct direction. This 
can lead to negative temperatures for the anisotropic thermal diffusion equation. In a previous paper we 
proposed a monotonicity-preserving explicit method which uses limiters (analogous to those used in the 
solution of hyperbolic equations) to interpolate the temperature gradients at cell faces. However, being 
explicit, this method was limited by a restrictive Courant-Friedrichs-Lewy (CFL) stability timestep. Here 
we propose a fast, conservative, directionally-split, semi-implicit method which is second order accurate in 
space, is stable for large timesteps, and is easy to implement in parallel. Although not strictly monotonicity- 
preserving, our method gives only small amplitude temperature oscillations at large temperature gradients, 
and the oscillations are damped in time. With numerical experiments we show that our semi-implicit method 
can achieve large speed-ups compared to the explicit method, without seriously violating the monotonicity 
constraint. This method can also be applied to isotropic diffusion, both on regular and distorted meshes. 
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1. Introduction 

Anisotropic diffusion equation arises frequently in diverse applications: microscopic transport in magne- 



tized plasmas (2|; image processing 17]; diffusion-tensor magnetic resonance imaging thermal properties 



of crystals [8|; transport in geological formations |5|, etc. In [22| we showed that simple finite-differencing 
of the anisotropic diffusion equation resulted in unphysical numerical heat fluxes, which lead to negative 
temperatures at large temperature gradients. Negative temperatures, in addition to being unphysical, result 
in an imaginary sound speed and associated numerical instabilities. 

Anisotropic diffusion equation satisfies important mathematical properties such as monotonicity preser- 
vation along the direction of diffusion, so that no new extrema are created in that direction and any existing 
extrema are not accentuated, i.e., maxima must drop or be unchanged, and minima must increase or be un- 



changed (e.g. see 11|.|14| and references therein). For simplicity we will refer to this as an extrema reducing 
property. In two or three dimensions the anisotropic diffusion equation of the form we are considering here 
can always be transformed to the form dT(a, /?, t)/dt = d/da[D(a, (3)dT/da\, where a is a coordinate that 
varies along field lines and (3 are field-line label coordinates that are constant along a magnetic field line, so 
the monotonicity properties of 1-D diffusion should be preserved. 



In [22| we proved that the use of slope limiters (e.g., see [12] ) to interpolate temperature gradients at 
cell faces guarantees that temperature extrema are not accentuated, as required physically. The extrema- 
reducing property ensures that the temperature is positive for a CFL stable timestep. Since the temperature 
is positive, numerical instabilities that plague simple finite differencing of anisotropic diffusion (because of 
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an imaginary sound speed!), do not arise with the use of limiters. Because of this desirable property our 
method has been used in astrophysical magnetohydrodynamic (MHD) simulations where thermal conduction 
is anisotropic, and large temperature gradients can arise [3, HH, 0|. In addition to thermal conduction, 
limiters have proved useful for anisotropic viscosity with large gradients [§| , and may help in various problems 
with anisotropic transport; e.g., cosmic ray streaming (23| . There is another explicit method that has been 
applied to astrophysical problems and avoids the problem of negative temperature in presence of large 
temperature gradients, though it is non-conservative and somewhat more complex We also discuss 

some other recent work below. 

A key limitation of all explicit methods is that the timestep is limited by the usual CFL condition, 
At < Ax 2 /2x\\, where Ax is the grid spacing and \\\ is the anisotropic diffusion coefficient. For some 
applications, this timestep constraint is rather severe and the conduction timestep can be much smaller 
than the MHD CFL time limit. In such cases an implicit method, where there is no stability limit on the 
diffusion timestep, is desirable. Although it is straightforward to difference a linear anisotropic diffusion 
equation implicitly, the resulting scheme is still not monotonicity-preserving. E.g., see Table 1 in 0], which 
shows that temperature oscillations remain till late times even with an implicit method, just as with the 
explicit schemes without limiters. One can try to solve the anisotropic diffusion equation fully implicitly, 
using limiters to prevent temperature oscillations, but the nonlinearities in the limiters will require a careful 
iterative treatment (some studies have found that these kinds of nonlinear limiters make iterative solvers 
more difficult). 

We have experimented with a Jacobian-free nonlinear iterative implicit method (a two-stage Richardson 
iteration extension of the LGMRES(1,I) version of Loose GMRES 2], a variant of the Generalized Minimal 
RESidual methodpH with restarting) but found that it requires a fairly large number of iterations per time 
step, because the Jacobian matrix is not strongly diagonally dominant for large timesteps and has a large 
condition number^] The fast method proposed here might be able to serve as an effective preconditioner to 
further accelerate an unsplit iterative method. 

We have also experimented with an explicit method which is stable for timesteps longer than the CFL 
limit, and where the internal iteration time-steps are chosen based on the properties of Chebyshev polyno- 
mials [lj. We were not able to obtain a speed-up of more than ~10 compared to the CFL-limited scheme, 
irrespective of resolution, for any of the parameters that we varied. Moreover, the parameters for which 
maximum speed-up is obtained, without becoming numerically unstable, are difficult to choose (this is true 
even for isotropic diffusion!). 

Here we present a conservative, directionally-split, semi-implicit method which is numerically stable for 
any choice of timestep, and is easy to implement. The method is based on directional splitting where the heat 
fluxes in each direction are updated sequentially. The heat flux in each direction, e.g., q x — —X\\bx(p ■ V)T 
(see Eq. [5]), consists of two terms: —x\\b x dT/dx, the 'normal' term where temperature gradient need not 
be interpolated, and —x\\b x b y dT/dy, the 'transverse' term which involves temperature interpolation with 
limiters. In our directionally-split method the 'normal' terms are treated implicitly, and the transverse 
terms are treated explicitly. The directional splitting of 'normal' implicit terms results in a tridiagonal 
matrix which can be solved very quickly. The explicit treatment of 'transverse' terms with limiters ensures 
that extrema are not accentuated. The resulting scheme, while not strictly monotonicity-preserving for large 
timesteps, results in only small amplitude temperature oscillations which are damped in time. Speed-up of 
order 100-1000, compared to the explicit scheme, is easily achieved for our test problems. 

There has been some interesting recent work on another approach to the problem of preserving positivity 
in presence of anisotropic diffusion tensors (or diffusion on distorted meshes), based on expressing the flux 
at cell faces in terms of the advected quantity at the cell centers on either side of the face (a "two-point flux 
expression"), but where the coefficients of this flux depend on the transverse gradients and so is nonlinear 
(for example see 13, 25|). Future work could compare the nonlinear limiters and implicit solvers that are 



2 Even for a linear 2-D Poisson problem on an TV X TV grid, Loose GMRES or conjugate-gradient methods require O(N) 
iterations by themselves, or 0{N 1 ^ 2 ) iterations if combined with a sufficiently good preconditioner like Modified ILU lldH . On 
the other hand, the splitting method employed in this paper uses tri-dagonal implicit solvers that are equivalent to only O(l) 
iterations and so are quite fast by comparison. 
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used in our algorithm with these other algorithms on the types of test problems considered here and in 22 1. 

The paper is organized as follows: Section 2 presents the method in detail and shows that the scheme 
is linearly stable for large timesteps. Section 3 presents results from three test problems which show the 
practical utility of our method. We conclude and discuss applications of our method in Section 4. 

2. The Method 

The anisotropic diffusion equation in its simplest form is given by 

dT 

at - -*'« (1) 

q = - X \\b(b ■ V)T = -X||W||T, (2) 

where T is the temperature, q is the heat flux along magnetic field lines, x\\ is the thermal diffusion coeffi- 
cient (with dimensions L 2 T _1 ), and b is the magnetic field unit vector. We will assume b as a given function 
of space and time, and \\\ as being a constant for simplicity (in [22| we showed that a harmonic average 
should be used for interpolation of the diffusivity for numerical stability). We assume two-dimensions with 
a uniform Cartesian grid and a constant diffusion coefficient. The generalization to a nonlinear diffusivity, 
general coordinate system, and three dimensions is straightforward. 

The semi-implicit method is obtained by directional splitting of the heat flux updates in each direction, 
given by 

jij i,j _ ^2 i,j _ ^2 i,j 



At Ax 2 x,i-l/2,j Ax 2 

, b x l+1 /2jby }i+1 /2,j -j^n fra:,i-l/2,jfti/,i-l/2j -fr7p" 

i,3 ~ i-J _ j2 ~ i,3 l2 i,j _ 

0„ ; ;xi ft : T\ 0„, 



XII At I/.W+V2 Ay 2 Ay2 

by,i,j+l/2b x ,ij+l/2 1 ^* ^J-l/2jgjj-l/2 -r7p* 

where magnetic field unit vectors are interpolated at the appropriate cell faces (simple averaging is fine for 
magnetic field unit vectors), and 

AT t+1 /2j = L (Tj+ij+i - Ti + ij,T i+ ij — Ti + ij^i,Tij + i — Tij,T hJ — Jjj-i) , (5) 
A^i,j+i/2 = L (Tj + ij + i - T hJ+ i,T hJ+ i — Ti_ij + i,Tj +1) j — Tij,T hJ — T^ij) , (6) 

are the temperature differences centered at appropriate faces, and L stands for a limiter. See [22| for a 



discussion of limiters in this context; here we will use the slope limiter of van Leer [27], If 21 ]; L(a, b, c, d) — 
L(L(a, 6), L(c,d)) is symmetric in its arguments, where 

L(a,b) = , if ab > 0, 

a + o 

= 0, otherwise. (7) 

We have experimented with other slope limiters (e.g., minmod and monotonized central [MC] limiters). We 
observed that monotonicity properties are much better with diffusive limiters (such as minmod and van-Leer) 
as compared to sharper limiters such as the MC limiter (see [12| for properties of different limiters) for our 
semi-implicit scheme with large timesteps. However, more diffusive limiters result in larger perpendicular 
diffusion. 
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Defining the components of the diffusion operator on the right hand side of Eq. ((T|) as 

(there is no implied summation on the right hand side of this definition), then the method in Eqs. ([3]) & 
dl} can be expressed as 

T n+i = (1 + AWyy yi^ _ AW yx )(l + AWxx)- 1 ^ - AW xy )T n . (9) 

Our formulation treats the 'normal' temperature derivative terms, which are guaranteed to result in 
a heat flux in the correct direction, implicitly. The 'transverse' terms are treated explicitly, and employ 
limiters that ensure that the temperature extrema are not accentuated. Directional splitting results in a 
quickly solvable tridiagonal matrix for each directional update. Instead of updating q x followed by q y we 
can also update the heat fluxes in the reverse order. Results do not depend substantially on the order of 
updates. Since our split operators are individually only first order accurate in time, Strang-splitting (e.g., 
see [lij p will not improve the accuracy of our method (and high accuracy is not a priority for components 
of the solution that are strongly damped anyway). Moreover, Strang-splitting applied to Eqs. ([3]) and (j4]) 
is numerically unstable for large timesteps. 



(8) 



2.1. Linear Stability Analysis 

One can perform the von Neumann linear stability analysis on Eqs. Q and Let us assume a single 
temperature mode T(x,y,t) = Tor(t)e~ l ( k:cX+kvV \ where r(t) is the amplification factor in time, and k x , k y 
are the wavenumbers in the x— and y— directions. The amplification factor can be written as r — r±r2, 
where ri and r<i are amplification factors for the substages in Eqs. and respectively. On substituting 
the discretized temperature eigenmode in Eq. (|3]), and using trigonometric identities, one obtains 

1 ~ sm(k x Ax) sm(k y Ay) 

r i = — X7 ; (10) 

1 + 4^6 2 sin 2 ^ Ax/2) 

where, for simplicity, we have assumed that b x , b y are constant in space. We use arithmetic averaging instead 
of the nonlinear limited averaging for the transverse temperature gradient. Similarly, for update in the y— 
direction, 

1- js^b x b y sm(k x Ax)sm(k v Ay) 

l + ^blsm 2 (k y Ay/2) ' " 
Now the amplification factor after a full timestep r = r\ri is 

[1 - A x A y cos{k x Ax/2) cos(k y Ay/2)} 2 



(1 + A*)(l + Al) 



(12) 



where A x = 2yfjqAlj^ s'm(k x Ax/2) and A y = 2^x\\At^ sm(k y Ay/2). From Eq. ((HJ) we get 

r< (1 + IA^I) 2 _1 + AIAI + 2\A x A y \ 



A2 A2 i A2 I A2 1 i A 2 A 2 i A 2 i A 2 ' 



which is guaranteed to be < 1 since 2,\A x Ay | < A 2 X + A 2 y for all real A x , A y . Thus our scheme (Eqs. 

SI is unconditionally stable like the usual implicit methods (see Fig. [IJ. The unconditional stability of our 

scheme also holds for the case of a general symmetric, positive diffusion tensor 



D = 



4 







X 



Figure 1 : Contour plots of the amplification factor \r\ = | r'i r2 1 (see Eqs . 1101 1116 for b x = b y = 1/ \/2 and different Courant factors 
(ncfl, see Eq. 1161 . Also shown is the analytic amplification factor (e k l\ Ax ncfl/4^ assumm g ^ x _ wnere fo,, = k x b x + k y b y ) 
for ncfi=10; there is no damping along k x = — k y because temperature gradient (which is along k) is perpendicular to b. 



with non-negative eigenvalues (i.e., d xx d yy > d xy ). Although large timesteps are stable, we cannot use very 
large timesteps because of loss of accuracy. 

Fig. [1] shows the amplification factor (|r| = | ri T2 1 ) for a timestep longer than the Courant time step 
(characterized by ncfl; see Eq. [T6|) for a fixed magnetic field unit vector b x = b y = \/\[2. As expected, 
the amplification factor is close to the analytic expectation (also shown in the same figure) for smaller ncfl. 
The damping rate is not sufficiently large (compared to the analytic solution) for unit vectors corresponding 
to the first and the third quadrant in /c-space. The maximum growth rate for the points in the first and 
third quadrants occurs at small scales (compared to the box size) , but this scale and the amplification factor 
becomes larger for a larger ncfl. Similar low damping rate arises for the Crank-Nicolson method for isotropic 
diffusion for large time steps. The slow damping of modes parallel to b at small scales can be avoided if we 
modify our scheme to a 4-step scheme, with the implicit terms in Eqs. & (H|) only applied for At/2 and 
two extra fully-implicit steps applied for At/2. Schematically this 4-step scheme is given by (see Eq. [5] for 
notation) 

T n+i = (1 + AWyy / 2 )- 2 (l - AW yx ){\ + AtV xx /2Y 2 {\ - AtV xy )T n . (13) 

However, for the test problems discussed in <J3]our 2-step method behaves satisfactorily even for ncfl as large 
as 1000. And moreover, perpendicular diffusion, non-monotonicity, and computational cost are slightly 
worse for the 4-step method (one can try different orderings of the operators in Eq. [T3] but this, and its 
analog where x and y are interchanged, worked best for our test problems). Thus we do not discuss the 
4-step method in detail. Another feature in Fig. [TJ which is seen for all ncfl, is that the smallest scale modes 
in the direction perpendicular to field lines are damped; this corresponds to cross-field diffusion at small 
scales. 

Expanding the amplification factor (r; Eqs. IT01 fc [TTj) in the limit k x Ax, k y Ay <C 1 (i.e., large length 
scales), one gets 

(1 - X\\k x k y b x b y At) 2 

r (i + X |i*2&2At)(i + xii*g&2Af)' 1 > 



which should be compared to the analytic amplification factor r a = e~ fc n X||A *, where fcii = k x b x + k y b y . In the 
limit x\\k 2 At <C 1, the difference between the numerical and analytical amplification factors is C([x|| k 2 At] 2 ); 
i.e., this method is first order in time. 

It is instructive to do a similar stability analysis in three dimensions, with n, r2, ^3 as the amplification 
factor in each directional update. In this case, 

1 ~~ SS^^A sin(k x Ax) siri(k y Ay) - j^-b x b z sin(k x Ax) sin(fc z Az) 
l + ^b 2 x sin 2 (k x Ax/2) ' 

and analogous expressions are obtained for other directions. It is easy to show that the absolute value of 
the amplification factor \r\ — {ri^r^] is not guaranteed to be < 1 for a large At. Thus, our scheme is not 
unconditionally stable in three dimensions. By numerically evaluating \r\ for different parameter^ we have 
verified that the stability condition for our scheme in three dimensions is (assuming Ax = Ay = Az) 

The corresponding stability limit for the explicit scheme in three dimensions is X||A£/A:r 2 < 0.444. Thus, 
our scheme can attain a speed-up of « 18 relative to the explicit method. Although our 2-step scheme (Eq. 
[9]) is not unconditionally stable in three dimensions, we have numerically verified that the 4-step method 
(Eq. II 3[) is unconditionally stable in three dimensions. 

We have also experimented with a variant of the alternate direction implicit (ADI) schemes for Eq. ([1]), 
inspired by their application to isotropic diffusion [lj| . Specifically, we tried (see Eq. [5] for notation) 

T n+i = (1 + AWyy / 2 )- 1 {l - At[D yx + V xy + V xx ]/2)(1 + AW xx /2)- 1 {\ - At[V xy + V yx + V yy ]/2)T n . 

However, presence of the transverse terms (d 2 T / dxdy) makes the scheme unstable for timesteps larger than 
a few times the CFL timestep, unlike in the case of isotropic diffusion where it is unconditionally stable. 
Even for isotropic diffusion, ADI does not give strong damping in the large time step limit, i.e., it is A-stablc 
but not L-stable. Therefore, we do not consider ADI further. 

Notice that our fully implicit scheme is only first order accurate in time, but for dissipative processes this 
is often adequate, since one is most interested in well-resolved components of the solution which are only 
weakly damped in a single time step. Another variant of Eqs. Q and (j4]), where the explicit 'transverse' 
{d 2 T / dxdy) terms are symmetrized with respect to the x— and y— updates (we were trying this to see if 
this scheme has better monotonicity properties as compared to our method), is numerically unstable for 
large Ats even though the linear stability analysis predicts an unconditional stability. Thus, linear stability 
is only a necessary (and not sufficient) condition for numerical stability, especially since the limiters are 
nonlinear. 



3. Numerical Tests 

In this section we describe various tests for our semi- implicit scheme (Eqs. [3] and 2]). 

3.1. Diffusion in a ring 

The ring diffusion test (see 0, H3|) involves the diffusion of a hot patch in fixed circular magnetic field 
lines. This is a crucial test to check monotonicity properties of the anisotropic diffusion scheme because 
field lines make all possible angles with respect to the Cartesian grid. At late times (a few diffusion times 



3 We calculate the maximum of |r| = |rir2T3| on a grid with resolution upto 320 3 in (k x Ax, k y Ay, k z Az), for a given 
{X\\Atb%/Ax 2 ,x\\ Atbl/ Ay 2 ,x\\Atb 2 z / Az 2 ). Then we try to find the maximum of X||At/Aa; 2 (assuming Ax = Ay = Az) for 
which \r\ < 1. 
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Figure 2: Temperature contour plots at t = 20 for the ring diffusion test problem using a 512 X 512 grid. Our semi-implicit 
method (Eqs. [3]and[4]| is used with different CFL numbers (ncfi). Also shown is the temperature plot with the fully explicit 
method (using the van Leer limiter) for comparison. 



across the ring), the temperature is expected to be uniform along each magnetic field line in the ring. The 
computational domain is a [— 1 , 1] x [— 1 , 1] Cartesian box. The initial temperature distribution is 

11 13 

T = 10 if .5 < r < 0.7 and — tt < 6 < —tt, 

= 0.1 otherwise, (15) 

where r = x 2 + y 2 and tan0 = y/x, and b x = —yj\/x 2 + y 2 , b y = xj 1 x 2 + y 2 . Reflective boundary con- 
dition (dT/dx — at boundaries in the x— direction, dT/dy = at y— boundaries) is used for temperature; 
magnetic field and conduction vanishes outside r = 1. The parallel conduction coefficient x\\ = 0.01; there 
is no explicit perpendicular diffusion. 

We quantify the timestep by ncfl, where 

At = ncflAa; 2 /4x||. (16) 

Since our scheme is unconditionally stable in two dimensions, we experiment with the CFL number (ncfl) 
to quantify the speed-up relative to the explicit method that we can obtain, without degrading the solution. 
Fig. [5] shows the temperature contour plots at t = 20 using different ncfl for a 512 x 512 box. As expected 
from section [2TT1 our scheme is numerically stable even for ncfl^ 1. However, the solution deteriorates for an 
extremely large ncfl; e.g., the temperature profile for ncfl=10000 looks quite different from rest of the others 
as there is considerable numerical diffusion out of the circular ring. This figure shows that large speed-ups 
(~ 1000 in the case of Fig. [2]) are possible as compared to the explicit method. Moreover, temperature 
oscillations at cxtrema are not as severe as schemes without limiters. 

Extrema are not accentuated for our semi-implicit method because the transverse term, responsible for 
non-monotonicity, vanishes at temperature extrema because of limited averaging. However, it is not guar- 
anteed that the temperature will be bound by the initial temperature extrema; this is because temperature 
oscillations may arise at non-extremal locations for a large CFL factor (ncfl). These newly created extrema 
will not be accentuated though. We never encountered temperature oscillations for the fully-explicit method 
with limiters which uses a CFL-stable timestep. 

Fig. [3] shows the temperature profiles at t = 20 for different grid resolutions but with a fixed ncfl=1000. 
The figure shows that the temperature profiles are very similar for n > 512. The maximum speedup relative 
to the explicit method (i.e., maximum value of ncfl), without seriously affecting the solution, is achieved 
for the highest resolution simulations (where significant speed-up is in- fact desired), as seen from Figs. [2] 
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Figure 3: Temperature contour plots at t = 20 for different grid resolutions but a fixed CFL number (ncfl=1000), using our 
semi-implicit method. 

& [3] While ncfl=1000 can be used for n = 512, the solution for a lower grid resolution with ncfl=1000 is 
quite diffusive in the perpendicular direction; also parallel diffusion seems to be suppressed (as is the case 
for ncfl=10000 and n = 512 in Fig. EJ. 

As mentioned earlier, our semi-implicit scheme is not guaranteed to be monotonicity-preserving. How- 
ever, the oscillations originating at large temperature gradients are of small amplitude and are damped 
away quickly in time (see Fig. 3]). In contrast, temperature oscillations for the explicit methods without 
limiters are large and persist till late times (see Fig. 7 in [22| )■ Of course, the amplitude of temperature 
oscillations is proportional to the ratio of maximum to minimum temperature at the discontinuity, but Fig. 
[4] shows that the temperature is still maintained positive for ncfl as large as 10000! The temperature ratio of 
100 (as in our test problem) is similar to the temperature range found in practical applications such as the 
transition of the solar chromosphere at 10 4 K to the coronal temperature of 10 6 K. Notice that the minimum 
temperature respects monotonicity constraint for the first timestep because all points in the initial condition 
are extrema and the transverse term vanishes. The temperature oscillations start from non-extremal points 
at later times. 

Fig. [5]shows the minimum temperature over the domain as a function of time for different grid resolutions 
and a fixed ncfl=1000. While Fig. 2] shows a clear trend of increased non-monotonicity as ncfl is increased, 
there is no systematic variation in the magnitude of temperature oscillations with the grid resolution for a 
fixed ncfl. This is expected because the factor x||At/Ax 2 is of the same order (~ncfl) for all resolutions. 
However, Fig. [3] clearly shows that for a fixed ncfl a more accurate temperature profile is obtained for a 
higher grid resolution. 

For a realistic problem ncfl should be < (I /Ax) 2 , where I is the scale on which we want temperature to 
be calculated accurately. Notice, that this factor increases with the grid resolution for a fixed I, so higher 
resolution runs will more accurate for a fixed ncfl (see Fig. [3]). Another constraint on ncfl comes from the 
positivity requirement. The magnitude of non-monotonicity is roughly independent of the grid resolution for 
a fixed ncfl (see Fig. [5]) , but depends on the ratio of maximum to minimum temperature at the discontinuity. 
E.g., for our test problem the initial temperature ratio is 100 and the maximum relative non-monotonicity 
(defined as {T m ; njU — min[T]}/T , m i njU , where T m i njU is the initial minimum temperature and min[T] is the 
minimum temperature of all times) for ncfl=1000 (see Fig. [5]) is s» (.1 — .08)/. 1 = 0.2. We have numerically 
verified that the relative non-monotonicity scales with the maximum to minimum temperature ratio (and 
of course non-monotonicity is larger for a larger ncfl). As mentioned before, non-monotonicity is worse for 
steeper limiters such as the MC limiter, but is less severe for diffusive limiters such as minmod. 
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Figure 4: Minimum temperature over the computational domain as a function of time (all timesteps are plotted) for our 
semi-implicit method using different Courant factors (ncfi). As in Fig. [2] grid resolution is 512 X 512. For comparison, the 
minimum temperature for the explicit method without limitcrs and with ncfl=l is -0.41. 




Figure 5: Minimum temperature over the computational domain as a function of time (all timesteps are plotted) for different 
resolutions using our semi-implicit method. As in Fig. [3] the Courant factor (ncfl) is fixed to be 1000. 
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Figure 6: The ratio of perpendicular (numerical) to parallel diffusivity as a function of grid size (triangles) for the smooth test 
problem in [2(|. The Courant factor (ncfl) is fixed to be 1000. Dotted line shows a second order convergence. 

3.2. Convergence & measuring xi, n um 

We perform a test problem with a smooth solution, described in 26] , to measure perpendicular numerical 
diffusion as a function of grid resolution and the Courant factor (ncfl). A two-dimensional Cartesian box 
([-0.5,0.5] x [-0.5,0.5]) is initialized with a zero temperature. Temperature is fixed to be zero at the domain 
boundaries at all times. We solve the anisotropic diffusion equation (Eq. [T|) with a source term 



where Q = 2ir 2 cos ttx cos iry. The fixed magnetic field is generated by a flux function <f> oc cos7rxcos7ry, so 
that the magnetic field unit vectors are along the contours of constant Q. And since temperature is driven 
by the source term, temperature is always constant along field lines. If there is no diffusion across field 
lines, temperature should rise with time. However, because of finite numerical diffusion in the perpendicular 
direction, it reaches a steady state. The steady state solution for the temperature, if we assume a finite 
perpendicular diffusivity \±> is T = Xj 1 cos 7rx cos 717/, independent of x\\- We use the asymptotic value (in 
time) of the maximum temperature to calculate X-L,num = 1/T(0,0). This test is slightly modified from 
[26j in that we do not include an explicit perpendicular diffusivity; this is because for the problems of our 
interest perpendicular conduction is negligible. 

Fig. [6]shows the ratio of the perpendicular numerical diffusivity and the parallel diffusivity (x±,num/x||) 
as a function of grid resolution for a fixed ncfl=1000. Perpendicular diffusion, which scales with x\\> shows 
close to a second order convergence with the grid resolution, as expected. Numerical diffusion is not sensitive 
to the Courant factor (ncfl); e.g., we verified that X-Unum/X|| is roughly independent of ncfl up to ncfl=10000 
for n = 256. Close to second order convergence of perpendicular numerical diffusion (and the independence 
of ncfl) was also seen for the ring diffusion test problem in section [3Tl 

3.3. Thermal Instability 

We also tested our method for a realistic astrophysical application, namely thermal instability in the 
intracluster medium, the X-ray emitting hot plasma pervading the massive clusters of galaxies. For as- 
trophysical motivation and details about numerical set up see [24| . We perform two-dimensional MHD 
simulations with anisotropic thermal conduction, using a periodic Cartesian box (40 kpc x 40 kpc) with 



dT 
~dt 



(17) 
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ncn=At/Ai eX p ncfl=l logic, r 




x(kpc) x(kpc) 

Figure 7: Logio temperature (keV) at 0.95 Gyr (pb 10 cooling times in the initial state) using our semi-implicit method (left) 
and the explicit method (right). The timestep for the conduction module in the semi-implicit scheme is chosen to be equal 
to the CFL timestep for the rest of the code (At) so that ncfl=At/At e xp, where At cxp = Ax 2 /4x|| jmax , and XlLmax ' s the 
maximum thermal diffusivity over the whole box. The conduction module is subcycled PS Ar,/At exp times for the explicit 



1024 grid points in each direction. The initial temperature is 0.78 keV and initial electron number density 
is 0.1 cm -3 . Identical pattern of small amplitude density /temperature perturbations are initialized (such 
that the pressure is uniform) to seed the thermal instability. The functional form of heating and cooling is 
such that cooling increases faster than heating for the cooler plasma, and vice versa for the hotter plasma. 
Thus, heating/cooling runs away and the plasma segregates into a two-phase medium. Net heating aver- 
aged over the whole box equals net cooling, so that the total thermal energy content of the box does not 
change with time. Thermal instability is in the isobaric limit; i.e., cooling time > sound crossing time over 
all scales. Thermal conduction is primarily along field lines aligned initially at 45° to the box; field lines 
roughly maintain their geometry even in the nonlinear stage. Small diffusion perpendicular to field lines is 
added for numerical convergence (see 24| for details) . Cold filaments aligned along the direction of the local 
magnetic field arise nonlinearly (see Fig. [7]) because thermal conduction along field lines suppresses growth 
of small scale modes. 

Fig. [7] shows the temperature in the nonlinear state of thermal instability obtained by treating thermal 
conduction using our semi-implicit scheme (left) and using the explicit scheme (right) with the van Leer 
limiter. The temperature plots are almost identical, establishing the practical utility of our method. The 
conduction timestep (At exp = Ax 2 /4x|| jlnax , where x\\, max is the maximum thermal diffusivity over the whole 
box and Ax = Ay is the grid size) in the initial state is 3 times the CFL timestep limit of rest of the code 
(At). As thermal instability becomes nonlinear and the hottest plasma becomes hotter, At oxp decreases 
rapidly relative to At, because of the sensitive dependence of conductivity on temperature (x\\ cx T 5 / 2 ; see 
Q). The temperature-dependent conductivity is interpolated at the faces using the current temperature 
(see for details about interpolation of conductivity). At 0.95 Gyr (the time corresponding to Fig. [7]) 
At/ At cxp ~ 100; thus the explicit method is subcycling the conduction module for 100 times, whereas the 
conduction module is applied only once for our semi-implicit scheme. Situation become worse with time 
because the hottest plasma in the box becomes hotter in time! Thus we are able to run much faster with 
our semi-implicit scheme, without affecting the solution and without violating temperature positivity. This 
example demonstrates the practical utility of our method. Also, recall that the stability limit for explicit 
\x 2 compared to Ax for the hyperbolic terms, so our scheme will be even more useful at 
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higher resolution and with mesh refinement. 



4. Conclusions 

We present a simple, directionally-split, conservative, semi-implicit method for anisotropic diffusion which 
is linearly stable for large timesteps. Directional splitting results in a tridiagonal matrix equation for each 
direction, which can be solved exactly and efficiently For problem on a N x N grid our scheme (Eqs. 
U]) requires two independent tridiagonal solves. In comparison, the fastest unsplit methods like LGMRES 
will require O(N) iterations to converge! Similarly, compared to the explicit method our scheme is ncfi (~ 
10-1000) times faster. Our method should be easily implemented in parallel using standard parallel linear 
algebra libraries like ScaLAPACK [(|. 

Although our method is not monotonicity-preserving for arbitrarily large timesteps, the temperature 
oscillations are of quite small amplitude and are damped with time. Using test problems we show that large 
speedups (up to ^100-1000 for our test problems) are achieved compared to the explicit method, without 
seriously violating the monotonicity constraint. A similar directional splitting may also prove effective for 
isotropic diffusion. Although ADI is numerically stable, fast, and second order accurate in time for isotropic 
diffusion, it does not give strong damping in the large timestep limit (just like the Crank Nicolson scheme; 
e.g., [H). 

We also tried unsplit methods, both fully implicit using limiters for the 'transverse' terms, and semi- 
implicit where only 'transverse' terms with limiters are treated explicitly. The limiters lead to nonlinearities 
that require some care with an iterative solver, and these unsplit methods result in a large sparse matrix 
equation which is much more expensive to solve than the tri-diagonal systems of the split method, even with 
an iterative solver like conjugate gradients or Loose GMRES. Although both of these unsplit methods with 
limiters appear to be monotonicity-preserving for arbitrary At, it takes many iterations to obtain a converged 
solution in both cases, and we generally find that the split algorithm is quite efficient by comparison. Our 
method might be able to serve as an effective preconditioner to further accelerate unsplit iterative methods, 
and works quite well as it is for our present purposes. 

Thermal conduction is primarily along the magnetic field direction for hot plasmas. For astrophysical 
plasmas large temperature gradients exist, and it is important for the numerical scheme implementing 
anisotropic conduction to yield positive temperatures in regimes of interest. Another practical requirement 
is that the scheme be fast so that the conduction timestep is not much smaller than the MHD timestep. Here 
we present a simple, directionally-split method that gets close enough to monotonicity that the temperature 
remains positive in presence of relatively large temperature gradients, and results in a large speedup. Our 
method should find applications in modeling of hot astrophysical plasmas with large temperature gradients 
(e.g., the multiphase interstellar/intracluster medium, transition from chromosphere to corona in the Sun). 
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